clear all;
close all;
delete(gcp('nocreate'));
%tic
%Define initial step for estimation algorithm;
global usual_delta;
usual_delta = 0.05;

%%%%%%%%%%;
%Folders;
%%%%%%%%%%;
% Current folder for Matlab codes
fp.matlab = pwd;
common_path = extractBefore(fp.matlab,"\matlab_dir");
fp.paper = sprintf('%s\\paper', common_path);
fp.boot_dir = sprintf('%s\\matlab_dir\\temp_boot', common_path);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We estimate the model at lambda=0.95;
fp.lambda = 0.95;

%model primitives;
fp.ages =9:12;
fp.periods=length(fp.ages);
fp.percentiles_kid = linspace(5,95,10);
fp.hc_points_kid = length(fp.percentiles_kid);

fp.Min_Time = 0.01;
fp.Max_Time= 1;
fp.inv_points=20;
fp.percentiles_inv= linspace(5,95,fp.inv_points);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%indexes for simulated data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.ind_data.child_id = 1;
fp.ind_data.child_age = 2;
fp.ind_data.n_sim = 3;
fp.ind_data.child_C = 4;
fp.ind_data.child_C_tp1 = 5;
fp.ind_data.inv = 6;
fp.ind_data.peers = 7;
fp.ind_data.peers_log = 8;
fp.ind_data.ParStyle = 9;
fp.ind_data.peers_tp1 = 10;
fp.ind_data.peers_log_tp1 = 11;
fp.ind_data.neighborhood = 12;
fp.ind_data.college = 13;

fp.ind_data_network.child_id=1;
fp.ind_data_network.child_age=2;
fp.ind_data_network.n_sim=3;
fp.ind_data_network.child_C=4;
fp.ind_data_network.n_friends=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loading Estimated Moments (and keep track of orders of moments with indexes);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

temp_mean_PS_between_age=importdata('mean_PS_between_class.txt');    
fp.mom_data.mean_PS_between_age = temp_mean_PS_between_age;

moments_temp=importdata('moments_educ.txt');    

fp.mom_data.mean_PS = moments_temp(1);
fp.mom_data.mean_PS_college = moments_temp(2);

fp.ind_mom.coef_PS.child_C = 1;
fp.ind_mom.coef_PS.peers   = 2;
fp.mom_data.coef_PS     = moments_temp(3:4);

fp.ind_mom.coef_skills_tp1.child_C    = 1;
fp.ind_mom.coef_skills_tp1.peers      = 2;
fp.ind_mom.coef_skills_tp1.PS         = 3;
fp.ind_mom.coef_skills_tp1.college    = 4;

fp.mom_data.coef_skills_tp1 = moments_temp(5:8); %Table 11 Column (1) 
fp.mom_data.coef_skills_tp1_ps0 = moments_temp(9:10);  
fp.mom_data.coef_skills_tp1_ps1 = moments_temp(11:12);

fp.ind_mom.coef_peers_tp1.child_C = 1;
fp.ind_mom.coef_peers_tp1.peers   = 2;
fp.ind_mom.coef_peers_tp1.PS      = 3;

fp.mom_data.coef_peers_tp1  = moments_temp(13:15); %Table 2 Column (2) 
fp.mom_data.coef_peers_tp1_ps0  = moments_temp(16:17); 
fp.mom_data.coef_peers_tp1_ps1  = moments_temp(18:19); 

fp.ind_mom.coef_inv.child_C = 1;
fp.ind_mom.coef_inv.peers   = 2;
fp.ind_mom.coef_inv.const   = 3;

fp.mom_data.coef_inv_PS0  = moments_temp(20:21); 
fp.mom_data.coef_inv_PS1  = moments_temp(22:23); 

fp.mom_data.mean_skills = moments_temp(24:27);

fp.mom_data.mean_inv    = moments_temp(28:29);

fp.mom_data.mean_nfriends = moments_temp(30);


fp.ind_mom.coef_PS_longit.child_C = 1;
fp.ind_mom.coef_PS_longit.peers   = 2;

fp.mom_data.PS_longit = moments_temp(31:32);
fp.mom_data.Inv_longit = moments_temp(33:34);

fp.mom_data.mean_PS_between = moments_temp(35:38);

fp.data_mom = moments_temp;

%Loading New Initial Conditions by Mother's Education;
initial_conditions = importdata('initial_conditions_educ.txt');    
fp.schools_distrib = initial_conditions;

%Loading Income information of neighborhoods (it affects preferences for parenting in the model);
mean_income_by_neighborhoods = importdata('family_income_neighborhoods.txt');
mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
fp.mean_income_by_neighborhoods = mean_income_by_neighborhoods;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%random mumber draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rng(10); %fix seed;
fp.n_rip_sim = 50;                          % # of simulation of the economy;
fp.tot_draws_mcintegration_network = 50 ;   % # of simulation of the network;

fp.init_draws{1} = norminv(rand(fp.schools_distrib(1,3),1,fp.n_rip_sim));
fp.init_draws{2} = norminv(rand(fp.schools_distrib(2,3),1,fp.n_rip_sim));
fp.init_draws{3} = norminv(rand(fp.schools_distrib(3,3),1,fp.n_rip_sim));
fp.init_draws{4} = norminv(rand(fp.schools_distrib(4,3),1,fp.n_rip_sim));

fp.peers_backward{1} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(1,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{2} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(2,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{3} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(3,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{4} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(4,3),fp.tot_draws_mcintegration_network, fp.periods-1);

%Here we double the shocks for peers formation in the backward problem
%because we double the state space compared to main problem by adding 
%mother's education (college vs. non-college)
for j = 1:1:4   
    fp.peers_backward{j} = [ fp.peers_backward{j} ; fp.peers_backward{j} ] ; 
end

fp.peers_forward{1} = rand(fp.schools_distrib(1,3),fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{2} = rand(fp.schools_distrib(2,3),fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{3} = rand(fp.schools_distrib(3,3),fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{4} = rand(fp.schools_distrib(4,3),fp.schools_distrib(4,3),fp.n_rip_sim,fp.periods-1);

fp.PS_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{4} = rand(fp.schools_distrib(4,3),fp.n_rip_sim,fp.periods-1);

fp.College_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim);
fp.College_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim);
fp.College_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim);
fp.College_forward{4} = rand(fp.schools_distrib(4,3),fp.n_rip_sim);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimate the model;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.do_estimation  = 0;
fp.do_counter_type = 0; 

if fp.do_estimation == 1

fp.do_counter = 0;            
                                                                                                                        
    load('param0')              
                  
    options_est = optimset;    
    options_est = optimset(options_est,'MaxFunEvals',3000);
    options_est = optimset(options_est,'MaxIter',3000);
    options_est = optimset(options_est,'TolX',1e+0);
    options_est = optimset(options_est,'TolFun',1e-2);
    options_est = optimset(options_est,'UseParallel',1);
    options_est = optimset(options_est,'Display','iter');
    param_est_educ = fminsearch_function(@obj_function, param0 , options_est,fp );
    save('param_est_educ','param_est_educ');
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%POST-ESTIMATION;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

do_post_estim = 1;

if do_post_estim==1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate:
%         1. Baseline Economy        
%         2. Print Sample-Fit
%         2. LateX Tables for Estimates;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('param_est_educ');   

fp.do_estimation =0;
fp.do_counter = 0;            
obj_eval = obj_function( param_est_educ, fp );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Create Table of new Estimates by Education;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

cd(fp.matlab) 
Estimates_tables_educ(obj_eval.structural_params,fp); 

end %end if do_post_estim = 1
